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Transient dynamics of an adiabatic NEMS 

M. Biggio\ F. Cavaliere^’^*, M. Storace\ and M. Sassetti^’^ 


This paper is focused on the transient dynamics of an adi¬ 
abatic nano-electromechanical system (NEMS), consisting 
of a nano-mechanical oscillator coupled to a quantum dot. 
By numerically solving the nonlinear stochastic differential 
equation governing the oscillator, the time evolution of the 
oscillator position, of the dot occupation number and of the 
current are studied. Different parameter settings are stud¬ 
ied where the system exhibits bi-stable, tri-stable or mono¬ 
stable behavior on a finite-time horizon. It is shown that, 
after a typically long transient, the system under investi¬ 
gation exhibits no hysteretic behavior and that a unique 


steady state is reached, independently of the Initial con¬ 
ditions. The transient dynamics is marked out by one or 
two well separated characteristic times, depending on the 
considered case (i.e., mono- or multi-stable). We evalu¬ 
ate these times for a dot on-resonance or off-resonance. 

It turns out that the characteristic time scales are long In 
comparison to the period of the uncoupled oscillator, par¬ 
ticularly at low bias, suggesting that the predicted transient 
dynamics may be observed in state-of-the-art experimen¬ 
tal setups. 


1 Introduction 

Nanoelectromechanical systems (NEMS) [1,2], consisting 
of a nanoscale oscillator coupled to the conduction elec¬ 
trons of a quantum dot, are intriguing systems from both 
an experimental and a theoretical point of view. Several 
different physical implementations of NEMS can be en¬ 
visioned. The oscillator is commonly fabricated through 
suspended nano-beams [3,4] or doubly-clamped carbon 
nanotubes (CNTs) [5-11], whereas the quantum dot can 
be either nano-fabricated in a semiconducting host or em¬ 
bedded in the nanotube itself. NEMS have attracted a con¬ 
siderable interest due to a wealth of possible applications, 
ranging from single-molecule mass spectrometers [12,13] 
to nanoscale gas sensors [14,15] and biosensors [16]. 
Basically, two markedly different regimes are possible for 
NEMS, according to the ratio between the bare oscillator 
frequency Ho and the average tunneling rate Eq of elec¬ 
trons flowing through the quantum dot. 

In the anti-adiabatic case (Qq » Eq), the oscillator motion 
is revealed by the presence of quantized vibronic side¬ 
bands in the conductance of the system [6,10,11,17,18]. 
One of the hallmarks of this regime is the Franck-Condon 
blockade, a peculiar low-bias current suppression occur¬ 
ring when electrons and vibrons are strongly coupled 


which has been theoretically predicted [19-23] and ex¬ 
perimentally observed [6,11]. Also peculiar features of 
the shot noise of conduction electrons have been pre¬ 
dicted [19,22,24,25] and the influence of quantum co¬ 
herences on the oscillator dynamics has been investi¬ 
gated [26,27], possibly leading to an effective intrinsic 
cooling even in the anti-adiabatic regime. 

In the opposite adiabatic case (Qq « Tq), which we in¬ 
vestigate in this paper, the bare oscillator dynamics is 
very slow in comparison to the electronic one. As a result, 
the oscillator behaves classically [28-31], with the ultra¬ 
fast motion of electrons giving rise to both an effective 
non-linear deterministic force acting on the oscillator and 
an effective non-linear damping [32-41]. Moreover, the 
fluctuating nature of electron transport in the quantum 
dot gives rise to a nonlinear stochastic force acting on 
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the oscillator [29,34,37]. This makes the system an ideal 
condensed-matter playground for the study of nonlinear 
oscillations [42-46] and determines a very rich dynami¬ 
cal scenario, subject of many theoretical studies. Owing 
to the nonlinear nature of the oscillator dynamics, phe¬ 
nomena such as hysteresis, multi-stability and switching 
have been studied [47-51]. Multistability in particular is 
a very intriguing feature envisioning NEMS as novel data 
storage elements. The key idea is to encode information 
in the different available states, to store data by preparing 
the system in one of the available ground states, employ¬ 
ing the electronic subsystem (e.g. suitably tuning a gate 
voltage) and to possibly read out the dot state at a later 
stage [32]. In this respect, a crucial parameter to asses the 
robustness of the information storage element is the aver¬ 
age dwell time of the states, i.e. the average time it takes 
for the system to switch between two different configura¬ 
tions. In this respect, the evolution towards a steady state 
may even be as detrimental since the stochastic jumps 
between different states renders the data storage useless. 
This fact motivates a thorough study of the transient dy¬ 
namics of an adiabatic NEMS, going beyond the large 
number of works which have on the other hand character¬ 
ized its steady state [28,29,31,32,37,41]. In addition, also 
the transient dynamics of a NEMS has been the subject of 
investigation [52-60]. Two closely related issues arise in 
this context, namely 

1. whether the system exhibits a unique steady state or 
not; 

2. how the steady state is approached throughout the 
transient dynamics. 

Both issues have been addressed theoretically and no 
unique answer has been provided so far. By employing 
time-dependent Hartree-Fock (TDHF) techniques when 
^ To [53], multi-configuration TDHF (MCTDHF) tech¬ 
niques when Qq « Tq [54] or even MCTDHF supple¬ 
mented by diagrammatic Monte Carlo techniques [55], 
the tendency towards a non-unique steady state, de¬ 
pending on the initial condition of the oscillator, has 
been pointed out. However, the above numerical analyses 
mainly focused on rather shorttime scales of just a few pe¬ 
riods of the uncoupled mechanical oscillator [55]. On the 
other hand, by employing the polaron tunneling approxi¬ 
mation for Do > Fo (thus, out of the adiabatic regime) [56] 
or in the presence of superconducting leads [58] no hys- 
teretic behavior has been reported. 

The results cited above evidence one of the main difficul¬ 
ties in studying the time evolution of NEMS: the presence 
of possibly very long time scales [55] makes it difficult 
to actually reach the steady state by means of numerical 
techniques. 


In this paper, we study the time-dependent evolution 
of the NEMS observables - oscillator position, dot level 
occupation and current - as a function of different sys¬ 
tem parameters. We focus on the adiabatic case and con¬ 
sider a strong coupling between the oscillator and the 
electrons, the most favorable condition to observe multi¬ 
stability [32], considering different possible experimental 
realizations and discussing their potential and limitations. 
Our task is to understand how the system evolves to a 
steady state and whether the latter is unique or not. De¬ 
scribing the system by means of the Anderson-Holstein 
model, the ultra-fast electronic dynamics is traced out, 
thus obtaining an effective nonlinear Langevin equation 
for the oscillator [32, 34] and a corresponding Fokker- 
Planck equation for the probability density of the oscilla¬ 
tor states [32]. By applying suitable numerical techniques, 
we solve the time-dependent Langevin equation and ob¬ 
tain the full dynamical evolution towards the steady state 
of the system observables with an arbitrary initial condi¬ 
tion. 

Our main findings are the following. 

1. Recasting the Fangevin equation in a Fokker-Planck 
form allows to prove that the system exhibits a unique 
steady state which does not depend on the initial con¬ 
ditions. 

2. When the system is in a multi-stable setting, the tran¬ 
sient dynamics is characterized by two or three time 
scales. The shortest ones are related to the dynami¬ 
cal trapping of the oscillator state around one of the 
equilibrium positions of the system. The longest one 
represents the characteristic relaxation time towards 
the (unique) steady state. 

3. When the system is in a mono-stable setting, the tran¬ 
sient dynamics is characterized by the relaxation time 
only. 

These characteristic times have been estimated by study¬ 
ing the statistical properties of the solutions of the 
Langevin equation and the spectral properties of the 
Fokker-Planck equation. Especially in the case of a multi¬ 
stable NEMS, they are found to be several orders of mag¬ 
nitude larger than the typical oscillator period I/Oq, sup¬ 
porting some of the claims made by other authors [55,56]. 
We remark that our numerical approach explores the 
system up to very long time scales, allowing to observe 
the fully developed steady state. In the multistable case, 
the shortest time scale is identified with the dwell time, 
which ultimately limits the performance of an hypotheti¬ 
cal NEMS data storage element. The characteristic time 
scales are found to be within the reach of state-of-the-art 
experimental investigations for many systems of interest. 
In particular, suspended CNTs seem ideal candidates in 
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view of their ability to reach the strong coupling adiabatic 
regime. 

The outline of the paper is as follows. In Sec. 2 the model of 
the system and the methods employed are presented. In 
Sec. 3 the results are presented and discussed, including 
possible realistic implementations of the NEMS under 
investigation. Some conclusions are drawn in Sec. 4. 


2 Model and Methods 


The NEMS is modeled as a single-level quantum dot lin¬ 
early coupled to a harmonic vibrational mode. The dot 
and oscillator Hamiltonians are [h - 1) 

Ha = ed^d, 


p2 mCll 

— -t- -X^ 

2m 2 


( 1 ) 


( 2 ) 


Here, e is the dot level position (which can be conveniently 
tuned via a capacitively coupled gate) and d is the dot 
fermionic annihilation operator. We assume here to deal 
with a dot characterized by a large charging energy, so 
that double occupancy can be neglected. Also, to keep the 
notation simple, a spinless model has been employed [32]. 
The extension of our results to the spinful case is straight¬ 
forward and does not lead to qualitatively different con¬ 
clusions. Position X and momentum P operators for the 
oscillator (with mass m) have been introduced, being Hq 
the bare oscillator frequency. The dot and the oscillator 
are coupled by the term 

Ha^o = ^XdU, (3) 


where A is the coupling force between electrons and vibra¬ 
tions. An example of the explicit form of A for the relevant 
case of a suspended CNT is provided in Sec. 3.4.2. The dot 
is also coupled to left and right contacts of free electrons 
described by 

Hi = ^ • ( 4 ) 

a=L,R k 

Ca,ic is the Eermi operator of lead a = L,R and k is the 
momentum of electrons therein. Contacts are kept at a 
given electrochemical potential ^a = Ro + sVa where e is 
the electron charge, Va a bias voltage applied to lead a 
and po the equilibrium chemical potential, assumed to 
be equal in both leads in the absence of bias. Henceforth 
we assume a symmetrical biasing with Vj, = -V12 and 
Vr = V/2. Dot and contacts are tunnel-coupled, with the 
tunneling Hamiltonian 

Ht= Y. YXadha,k + '^-c.. (5) 

a=L,R k 


The above Hamiltonians constitute the Anderson-Holstein 
model. This deceptively simple model produces a very 
rich and interesting physics and has been successfully em¬ 
ployed to describe NEMS in a vast range of regimes [19- 
21,29,32,34,49,53]. The dot-oscillator coupling sets the 
characteristic polaron energy and length scales 


Ep — 


A" 

2mkl^ 



( 6 ) 


whereas the tunnel coupling sets the average tunneling 
rate To = 27r|;|;|^ where, for simplicity, we assume symmet¬ 
ric tunnel barriers Xa = E,R)- 

In the adiabatic assumption (Oq « Tq), the vibrational 
dynamics follows adiabatically the ultra-fast motion of 
electrons, which act as an effective force and a dissipative 
bath on the oscillator [29,32,34,37,61]. Standard path- 
integral techniques [29,34,61,62] allow to trace out the 
degrees of freedom of electrons in both the contacts and 
the quantum dot, thus leading to an effective Langevin 
equation of motion for the oscillator 

jc -t A(x}x- Fix) = VD{x]^M , (7) 


where x = XI£p is the dimensionless oscillator posi¬ 
tion (regarded as a classical variable) and overdots imply 
derivatives w.r.t. t = Oq t, the dimensionless time variable. 
The dimensionless Eq. (7) is characterized by a non-linear 
and positive definite friction coefficient A{x), an effective 
potential U{x) with F{x) = -dxU[x] and a multiplicative 
noise term 

The nonlinear force term F (x) is due to the coupling of the 
oscillator to the dot and arises at the lowest (zero-th) or¬ 
der in D. 0 /T 0 (i.e., in the limit m — 00 ) [37]. The nonlinear 
damping and noise terms, on the other hand, represent 
the first-order contributions in Qo/ko to the oscillator dy¬ 
namics [37] (adiabatic corrections) and are related via the 
fluctuation-dissipation theorem [34]. In particular, the 
stochastic nature of the ultra-fast electron motion gives 
rise to random fluctuations of the dot occupation. These 
in turn are represented by the stochastic forcing term 
\/D{x]f{t), where ^(t) is a white Gaussian noise [29] with 
{fir)) - 0 and ((^(t)^ (t'))) = 5 ( t - t '). Here <<...>> denotes 
an average over the realizations of the stochastic process 

^T). 

In this work, we do not consider explicitly extrinsic damp¬ 
ing mechanisms due to possible external damping baths. 
Their expected impact on the results which we will present 
below are briefly addressed in Sec. 3.4.1. 

The dimensionless damping, force and noise terms are 
expressed by means of the forward and backward dot 
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Green’s functions on the Keldysh contour G+ [cj, x) as [32] 


— f da>' 
2jiJ 

G+ (a>', x)da,' G- (at', x), 

( 8 ) 

1 h; 

1 CM 

-t 

1 

^ dat' G+(a)',x), 

(9) 

— f da>' 
2jiJ 

G+(a>',x)G-(a)',x), 

( 10 ) 


where w = D.o/{2Ep). At small temperatures « Fq), 
the Green’s functions are [32] 


G±{cj,x') - ±ij ^ 




Xr {W-£-X)^+Y^ ’ 


( 11 ) 


The prohability density v; t) is obtained numerically 
by solving the SDE for a fairly large number of different 
realizations of the stochastic process ^v(t) by means of 
a highly optimized parallel algorithm. Solution traces for 
x{t) and v{t) are then sampled and a histogram is created 
for &>{x, v,r]. The process is iterated until convergence 
on the probability density is achieved. We remark that 
this procedure is not restricted to the asymptotic (t — oo) 
case. Indeed, as already anticipated, we are particularly 
interested into the transient evolution of the NEMS. 

An alternative approach to obtain &>{x, v, r) can be pur¬ 
sued, casting the Langevin equation (7) into a Fokker- 
Planck equation [63] ^(x, v, r) = [3^{x, v,r]], where the 

operator is given by [32] 


with 7 = ro/(2£p), £ = e/lEp, and d[x] the Heaviside step 
function. By substituting Eq. (11) into Eqns. (8-10) one 
obtains [32,37] 


A[x) = 

U{x) = 


2(1) 




(1-t A2 ) 
x(x-i-l) 7 


-y [Aa arctan (Aq 

2 i^ak.R 


( 


\ 


+ log 


1 + A^; 


(X) 


Dix) = — y 5 a 

a=L.R 


arctan (Aa) + ■ 


A„ 


1 + A^ 


( 12 ) 


(13) 

(14) 


(15) 


where 5 l = -tl, 5 ^; = -1 and 


Ai = 


U + 2[£- x) 

r 


Ar 


u-2{£-x) 

r 


(16) 

(17) 


with u = lelt7/(2£p). The shape of the potential energy 
U{x) (and then the oscillator dynamics) changes accord¬ 
ing to £ (related to the normalized gate potential) and u 
(the normalized voltage bias between the leads), which 
are chosen as bifurcation parameters. 


d 

5£ ]^] = - V— - F[x]^—+ A{x]—+D[x)—^. 

ox ov ov ov^ 

(19) 

Botb methods are employed here, since they bring com¬ 
plementary information on the system behavior. 


Once ^{x, v; r] is obtained, the expectation value of any 
given observable &{x, v] at time t can therefore be eval¬ 
uated as {©’(t)> = f dx f dv 0^[x, v;r)&[x, u). When these 
quantities do not depend on v but only on x, these ex¬ 
pression simplify to P{x\r) = f dv SP[x, v,r) and {©{r)) = 
f dx P(x; T)©ix). 

In this paper, we focus on three quantities, namely the 
position X, the dot current and its occupation number, 
given for fcg F « Tq by [37] 


i V- 

I[x) = — 2^ Sa arctan (Aa) 


2^ a=L,R 


( 20 ) 


n{x] = 


1 1 

- 1 - 

2 2ji 


y arctan (Aa). 

a=L,R 


( 21 ) 


The latter quantities depend parametrically on x in the 
adiabatic regime [32,37], while the time dependence of 
their averages is encoded in ^(x; t). 


3 Results 


Equation (7) is a stochastic differential equation (SDE). 
Therefore, given the v-th realization of the noise pro¬ 
cess fy(T), the solution Xy(T), subject to initial conditions 
Xy (0) = xo and Xy (0) = vo, fluctuates stochastically. A prob¬ 
ability density B^ix, v, r) for the oscillator to occupy the 
state (x, v) (with u = x) at normalized time r can be intro¬ 
duced as 

g^{x, v;t] - {{6 (x- Xy(T))5(t'- Xy(T))» . (18) 


3.1 Effective potential landscapes 

In the following, we will consider the strong-coupling 
regime ]32] 

7<1 (22) 

which, together with the adiabatic hypothesis a) « 7 to 
satisfy the adiabatic condition), restricts a)« 1. As we will 
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see, this is the most favorable regime to observe multi¬ 
stability and switching phenomena [32]. From the physi¬ 
cal point of view, the strong coupling regime corresponds 
to setting Ep as the largest energy scale. 

The equilibrium conditions for the noiseless system (Eq. 7 
with D{x) = 0) are v = Q and F{x) = -d^Uix) = 0. Then, 
the equilibrium points are the extrema ofU[x]. We start 


0.05 

U{x) 


0 

0.05 


U{x) 


0 

0.1 


U{x) 


- 0.1 



-1.3 


0.3 -1.3 


0.3 


Figure 1 Effective potential U{x) for different values of e and 
u. Left column: on-resonance case, e = -0.5. Right column: 
off-resonance case, e = -0.4. (a,d) u = 0.2; (b,e) u = 0.475; 
(c,f) u = 1.0. The red (blue) shade marks the regions where 
A{x] (D{x)) in Eqns. (12,14) is maximal. Other parameters: 
to =10“^, 7 =0.08. 


in shaping U{x) and giving rise to multistability, let us 
consider the most favorable case, that of a NEMS in res¬ 
onance e = -0.5 in the low bias regime u « 1. Here, the 
potential barrier AH = [/(-1/2)-H(0) = [/(-1/2)-[/(-I), 
separating the two degenerate minima, can be estimated 
as 



arctan|-j -i-ylog 


y 

\/i+r^ 


1 

8 ’ 


(23) 


One finds that for y < 1 one has AH > 0, while increasing 
7 above that threshold results in a progressive reduction 
of AH and the ultimate disappearance of bistability. 



discussing the shape of H(x), by varying the parameters 
u and e. The effective potential H(x) is shown in Fig. 1 
for an on-resonance (e = -0.5 [33], panels a-c) and off- 
resonance {£ = -0.4, panels d-f) dot. At low u values, two 
minima appear, around x = 0 and x= -1, for both the 
on-resonance (panel a) and the off-resonance (panel d) 
cases. Increasing u, the on-resonance case (panel b) is 
characterized by three minima - the two discussed above 
and a new one developing at x = -0.5. On the other hand, 
the off-resonance effective potential can have either two 
or three minima. Panel e displays a limit case with two 
minima and a third “ghost” minimum (see also Fig. 2b). 
For u = 1 or larger, only the minimum at x -0.5 survives 
for both the on-resonance and off-resonance cases (pan¬ 
els c and f). 

To illustrate the importance of the strong coupling regime 


Figure 2 Position of stable (solid line) and unstable (dashed 
line) equilibrium positions for the system as a function of the 
normalized bias u for (a) e = -0.5 and (b) e = -0.4. The green 
lines represent the values of u at which the plots in Fig. 1 have 
been made. 


Figure 2 shows the position of the minima (repre¬ 
sented by a solid line, and corresponding to stable equilib¬ 
ria for the noiseless NEMS) and maxima (represented by 
a dashed line, and corresponding to unstable equilibria) 
of H(x) as a function of u. For low values of u, in both 
the on-resonance (panel a) and off-resonance (panel b) 
cases, the presence of a bistable region is evident. As u 
increases, in panel a a subcritical pitchfork bifurcation 
occurs and the system enters a tri-stable region. Two fold 
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bifurcations at still higher u values mark the end of the tri¬ 
stable region and the beginning of a mono-stable region 
with only the minimum located at x = -0.5. Tuning the 
system off-resonance by changing e induces a symmetry 
breaking in the diagram described above. When moving 
far enough from resonance, the tri-stable region even dis¬ 
appears. Indeed, panel b shows the limit case (e = -0.4) 
where the tri-stable region disappears. 

There is a close relationship between the oscillator and 
the dot states. In particular, for given u and e, each mini¬ 
mum of U{x] corresponds to a different behavior of the 
quantum dot. The two minima at x =; 0 and x =! -1, oc¬ 
curring for low to intermediate u values, correspond to 
situations where the dot is essentially locked into either 
the empty [n = 0 ) or the occupied {n = 1 ) charge state re¬ 
spectively. To exemplify this, let us consider the low-bias 
case (u =2 0 ), where 


suggests that the system may reach a dynamical steady 
state characterized hy stochastic jumps between min¬ 
ima while the occupation prohahility v,r) attains 
a steady shape independent of the particular initial con¬ 
dition. This fact will be proven in the next section. 

3.2 Uniqueness of the steady state 

The adiabatic approximation considered here leads to the 
Fokker-Planck operator Eq. (19) dehnes a steady Fokker- 
Planck equation. For this equation, a unique steady state 
exists, for any initial condition [64]. To prove this fact, let 
us introduce the quantities 


V 

: a = 

Fix] - vAix] 

; d = 

' Dix] o' 





0 

0 



1 

1 

<2£\ 




niO] = 

-h 

— arctan 



(24) 

which 


2 

n 

7 1 



lows 


1 

1 

2 (e+l) 




ni-1] = 

-h 

— arctan 


(25) 



2 

71 

7 



if = - 

For £ > 

-1 and 7 « 1 

, one has n( 0 ) 0 and n(- 

1 ) = 1 . 



d 


u'- 


<J2 


2 dyf,dy. 




(28) 


This situation is the semiclassical counterpart of the well 
known Coulomb blockade regime [32]. 

The third stable minimum only arises for u> u*{£). For 
instance, for e = -0.5 one finds u*(-0.5) y/ZjIn. In 
general, as exemplified in Fig. 2(b), one finds that u* (e 7 ^ 
-0.5) > u* (-0.5). An approximation for the position of 
this minimum is 


where q,v- 1,2 and a summation over repeated indices 
is implied. In order for the Fokker-Planck equation to be 
steady and exhibit an unique solution, two conditions 
have to be met [64]. First of all, the matrix d must con¬ 
tain a positive-definite sub-matrix, which is obvious since 
Dix] > 0. The second condition is that the system of par¬ 
tial differential equations 




7 (2c -t 1) 

71 [u'^ + Y^) - 2j 


(26) = 0 

dx 


(29) 


valid for u>u* (e) and [c - 1 -1/2| < 7 . Roughly speaking, for 
u > 1, X* =2 -0.5. Additionally, it can be checked numeri¬ 
cally that n [x* (e)] =2 0.5 (with n [x* (-0.5)] = 0.5) and that 
for the same values the current exhibits a maximum. This 
situation is hence closely reminiscent of the sequential 
tunneling regime, where for symmetric tunnel barrier one 
expects half dot filling and maximal current through the 
system at hnite bias [32,33]. 

What described above holds for the noiseless system. The 
presence of the electronic noise induces crucial modifica¬ 
tions in the above simple picture. In particular, it induces 
jumps between the different minima of U{x) in the bi¬ 
stable and tri-stable cases, and stochastic fluctuations 
around the minimum of U{x] in the mono-stable case. 
The regions where D{x] is largest are marked with a blue 
shade in Fig. 1, whereas the regions where A(x) is maxi¬ 
mal are marked by red shades. 

Due to these jumps, the oscillator cannot remain indeh- 
nitely stuck into one of the potential wells. This already 


ds[x,v;r) dgix,v;T) 

^ -lF[x]-vA(x)]- ' 


(30) 


dr dv 

has the unique solution g(x, v,r) = const. [64]. This can 
be easily proven by further deriving Eq. (30) with respect 
to X and taking into account Eq. (30), which implies 


ds,{x,v,r] 

{d^F{x) - vd^Aix]] - = 0. 

dv 


(31) 


Since the quantity within square brackets is in general 
non-vanishing, one can conclude that di,g{x,v,r) = 0 , 
which in turn implies (owing to Eq. (30)) djg(x, v;r) = 0 
and shows that g(x, v; r) is constant. 

This proves that our Fokker-Planck is steady, from which 
stems that there is a unique steady state regardless the ini¬ 
tial condition of the system [64]. Therefore, any hysteretic 
behavior can be ruled out, in contrast to previous results 
obtained by means of mean field methods [53-55]. We 
want to stress here that the nature of this steady state is 
dynamical: indeed, when the steady state is reached one 
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has that0^(x, v, t) is independent of t but x{r] and v{tau) 
stochastically fluctuate due to the electronic noise term 
in the Langevin equation. 

Having proven that the steady state is unique, we still have 
to understand how the system approaches this unique 
steady state studying the full time evolution of the oscilla¬ 
tor and dot variables, which constitutes the main task of 
this paper. 


3.3 Convergence to the steady state 

For a given set of dot and oscillator parameters, we have 
solved Eq. (7) starting from different initial conditions for 
the oscillator and following the evolution of oscillator po¬ 
sition X, dot occupation n{x) and current I{x) (averaged 
over the noise realizations), until a steady state is reached. 


3.3.1 On-resonance case (£ = -0.5) 

Figure 3 shows plots of{x(T)), {/(t)> and(n(T)> for two dif¬ 
ferent initial conditions in the bi-stable region with bias 
u = 0.2. All curves converge in the very long time limit to 
a common steady state. We stress that, in accordance to 
what has been proven in Sec. 3.2, there is convergence to a 
unique steady state in all cases that we have investigated. 
The two initial conditions chosen here set the oscillator at 
rest near either of the minima of Ulx) for t = 0 . 

The dynamics shows two very distinct phases. Initially, the 
oscillator evolution is essentially/fozen within the poten¬ 
tial wells. Indeed, a close-up of this region, shown in the 
inset of panel a for (x(t]}, displays a damped oscillatory 
behavior with a quasi-period given by the bare oscilla¬ 
tor frequency. In this situation the dynamics is ruled by 
a competition between friction and noise (see Fig. 1(a)). 
Indeed, the damped oscillations are due to the weak tails 
of the nonlinear damping A(x). However, since the noise 
term is even smaller than A(x) near the minima of U[x), 
the oscillator can settle near the minima. While the os¬ 
cillator is locked near the potential minima, the dot is 
in a well-defined average state, empty or full, depending 
on the occupied well. This agrees with the discussion in 
Section 3.1 and shows how preparing the oscillator in a 
certain state correspondingly locks the dot charge. Only 
after a certain time, the small noise terms are able to in¬ 
duce transitions between the two minima. This marks the 
onset of a new system dynamics, dominated by stochastic 
jumps between the potential minima, induced by the noise. 
As soon as this transition sets in, the system begins to 
evolve towards the unique steady state, which is eventu¬ 
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Figure 3 (a) (x(t)>; (b) (/(t)>; (c) {n{r]) for u = 0.2, e = 
-0.5 and different initial conditions xq = 0, fo = 0 (red curve) 
and xo = -1, iiQ = 0 (blue curve). The inset in panel a shows a 
zoom of the short-time behavior of {x(t)> for the initial condition 
Xo = 0,1'o = 0. The orange (yellow) vertical line marks the 
characteristic dwell (relaxation) time scale (Xr) (see text for 
details). The steady state values are denoted by a dashed line. 
Other parameters: w = 10“^, j = 0.08. 


ally reached within a given time scale. Notice that the time 
scale needed to reach the steady state does not depend on 
the initial condition of the oscillator. In the steady state, 
the current is larger than in the transient evolution, due 
to the (partial) occupation of the region around x = -0.5, 
where the current reaches its maximum value. It has to 
be pointed out that the envelopes of (x(t)>, (/(t)), {n{r]) 
cannot be simply fitted by a simple exponential function 
(not shown), confirming the multiple time-scale dynam¬ 
ics of the system. 

We now dehne the two relevant time scales which have 
been described above and which characterize the tran¬ 
sient dynamics towards the steady state. We define dwell 
time Trf the average time that the system spends in the 
"frozen" configuration, and relaxation timexr the typical 
time scale over which the system reaches the steady state. 
Of course, < t^. As we will see later, may depend on 
the initial condition of the system and even collapse to 
zero in certain situations. On the other hand, is always 
different from zero and essentially independent of the 
initial conditions. 
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No multi-stability is hence found. However, in the case 



Figure 4 Relaxation time Xy (see text) as a function of u 
for e = -0.5 (red curve) and e = -0.4 (blue curve). The dots 
correspond to the cases shown in Figs. 3, 3, 7, 8, 9. Other 
parameters: t) = 10“^, y = 0.08. 


To determine a statistical analysis of several different 
solutions of the Langevin equation for different initial 
conditions and realizations of the noise, has been carried 
out for times up to t = 10®. For a given realization of the 



discussed above both and Xy are very long in compar¬ 
ison to the bare period of the uncoupled oscillator. This 
suggests that the relaxation towards a steady state can 
only be observed if an experiment (or a calculation) is per¬ 
formed up to very long times in comparison to the bare 
oscillator period. 

The two time scales are obtained analyzing the long-time 
behavior of the solutions x(t) of Eq. (7) and of the opera¬ 
tor ^ in Eq. (19). 

The relaxation time scale Xy is inferred by a numerical 
analysis of the Eokker-Planck equation. Indeed, the ex¬ 
istence of a unique steady state ensures the presence of 
only one zero eigenvalue, Cq = 0> of the operator S£. All the 
other eigenvalues have negative real parts. The eigenvalue 
with the smallest nonzero real part, e \, sets the longest 
time scale of the system. We have always found that e\ 
is real, and defined Tr = le^ ^ |. In order to determine this 
eigenvalue, has been discretized and the eigenvalue 
problem 5£v = Xv solved, using suitable algorithms for 
sparse matrices. Stability against the discretization of the 
operator has been checked. Figure 4 shows the value of Tr 
as a function of the bias u, in both the on-resonance and 
off-resonance cases. In both cases, Xy decreases as u is 
increased, and eventually saturates. Notice that the thresh¬ 
old value of u for this saturation roughly corresponds to 
the rightmost fold bifurcation(s) in Fig. 2, i.e. to the en¬ 
trance in the mono-stable region. 


Figure 5 (a) A typical solution trace Xy{x) (green curve), the 
running-averaged version (blue curve) and trigger-detection of 
the occupancy of oscillator minima (red curve), (b) Probability 
riofT) that the oscillator has spent a time x in the minimum of 
U(x) around x= 0. Parameters: u = 0.2, e = -0.5, w - 10~® 
and y = 0.08. 


noise ^v(t), the corresponding solution Xv(t) is smoothed 
by a running-average method. Subsequently, the regions 
where |Xv - x^^'' | < 5, |Xv - xj“^' | < 5 and |Xv - x|^* | < 5 are 
Identified (trigger detection), where x^^^ is the location of 
the minimum of the potential well corresponding to an 
average dot occupation k{k=Q, 1,1/2), and 5 is a thresh¬ 
old, set to 0.03 in our numerical analyses. This allows to 
determine the time x^J^ spent during the i -th visit of well 
k. From a statistical analysis of the times the prob¬ 
ability distribution flfclT) that the oscillator has spent a 
given time x in the fc-th well is obtained; the procedure 
is then iterated over several different realizations of the 
noise process. The procedure is exemplified in Fig. 5(a) 
and the resulting probability distribution for the well cen¬ 
tered at Xq"'' is shown in Fig. 5(b). Finally, a dwell time 
T® = f dxx IlfclT) is evaluated. In the on-resonance case, 
one finds = x^, whereas off-resonance two dif¬ 

ferent dwell times are obtained. 
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Figure 6 Reduced probability density ^{x;t) as a function 
of X for w = 0.2, £ = -0.5 and different values of t, given the 
initial condition xq = 0, vo = 0: t = 500 (red); r - 2-10^ (green); 
T = 5 ■ 10^ (cyan); t = 10® (blue). Other parameters: w - 10“®, 

j- 0.08. 


The orange and yellow vertical lines in Fig. 3 denote 
the estimated values for the dwell and relaxation times 
Td ~2 - 10® and Tr 4 ■ 10^, respectively. Figure 6 shows 
snapshots at different times of the reduced prohahility 
density of the oscillator with initial condition 

xq = 0,Vo = 0, one of the two cases shown in Fig. 3. For 
T < Trf the occupation of the well about x = -1 is clearly 
negligible, and a stable steady state is reached only for 
T > Tr (cyan and blue curves). Notice that the steady state 
density ^{x;t) is symmetric with respect to x = -0.5, 
as implied by the symmetry of U(x), Aix], and D{x) for 
£ = -0.5. 

Figure 5(a) also proves our statement about the dynamical 
nature of the steady state: even in the long time regime, 
the oscillator position stochastically jumps between the 
minima of the effective potential and only the probability 
distribution satisfied ^(x, v; r). 

The results in the tri-stable case are shown in Fig. 7, when 
the system is initialized near each of the three minima 
of U{x). The behavior when starting at rest from xq = 0 
and Xq = -1 (red and blue curves) is qualitatively simi¬ 
lar to that for the bi-stable case, with shorter and t^. 
The solution when starting near the third, middle, well 
is different (green curve). Indeed, both average position 
and dot occupation vary only slightly, whereas the current 
exhibits a marked decrease in time. No dwell time can be 
detected for this case, since, due to the noise term D{x], 



Figure 7 (a) (x(t)>; (b) {/(t)>; (c) (n(T)> for u = 0.475, 
£ = -0.5 and different initial conditions Xq = 0, Vq = 0 (fed 
curve), Xo = -1, fo = 0 (blue curve), and Xq = -0.45, Vq = 0 
(green curve). The orange (yellow) vertical line marks the char¬ 
acteristic dwell (relaxation) time scale (t^). The steady 
state values are denoted by a dashed line. Other parameters: 
oj = 10^®, 7= 0.08. 


the state quickly escapes from the central valley of U (x) 
towards one of the lateral minima, see Fig 1(b). 

The current traces are identical when the minima near 
x = 0orx = -1 are initially populated, similarly to the 
bi-stable case. The initial population of the central po¬ 
tential well (green curve) makes the current maximal, as 
discussed in Sec. 3.1. The current decreases as far as the 
lateral wells get occupied. We note that the decay time of 
(/(t)> roughly corresponds to Tr, which shows that for this 
kind of initial conditions the only time scale relevant to 
define the transient towards steady state is the relaxation 
time. 

For u > yjljln, only the minimum at x = -0.5 survives. 
Figure 8 shows the average position, current and dot oc¬ 
cupation versus t: in this case no dwell time can be de¬ 
fined as the dynamics is entirely restricted to the single 
central potential well. Indeed, all the solutions for the 
three initial conditions display essentially the same qual¬ 
itative behavior and the calculated value of Tr roughly 
matches the time when the transition towards the steady 
state occurs. Only the current trace for the initial condi- 
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Figure 8 (a) (b) (c) for m = 1, e = -0.5 

and different initial conditions jcq = 0, i/q = 0 (red curve), jcq = 
-1, 1^0 = 0 (blue curve) and xq = -0.45,1'o = 0 (green curve). 
The yellow vertical line marks the characteristic relaxation time 
Tr. The steady state values are denoted by a dashed line. 
Other parameters: w = 10^^, y = 0.08. 


tion jco = -0.5, Vo = 0 attains larger values at short times, 
since in this case the probahility distribution is initially 
concentrated near the region where current is maximal. 
The oscillations shown in panels (a) and (c) are due to 
the small number of points considered on the time axis, 
which does not allow to sample in full details the short- 
time oscillations already mentioned above. 


3.3.2 Off-resonance case (e = -0.4) 

Figure 9 shows the case u = 0.2 (see Fig. 1(d)). The most 
striking feature, as already anticipated in Sec. 3.1, is the 
presence of two different dwell times corresponding to 
the two different potential wells, with the most stable one 
having the largest t^, as expected. The average steady 
position decreases towards x = -1 and, correspondingly, 
the occupation of the dot tends to a value larger than 0.5, 
since in this case the steady state cannot be expected to 
be symmetrical with respect to x = -0.5. More intrigu¬ 
ing is the behavior of the current as a function of time, 
when the oscillator is initialized around x = 0, namely in 
the least favorable well. Here, the current exhibits a non- 


Figure 9 (a) (x(t)>; (b) {/(t)>; (c) {u(t)> for u = 0.2, e = 
-0.4 and different initial conditions xo = 0, yo = 0 (red curve), 
xo = -1, yo = 0 (blue curve) and xq = -0.45, uo - 0 (green 
curve). The orange (yellow) vertical line marks the characteris¬ 
tic dwell (relaxation) time scale (Tr). The steady state values 
are denoted by a dashed line. Other parameters: w - 10~^, 
7 = 0.08. 


monotonous behavior with a sharp increase leading to 
a maximum for r = r* ^ 5500, followed by a decrease 
towards the asymptotic value. In order to understand 
this behavior, in Fig. 10 we show the reduced probabil¬ 
ity density ^{x; t) for different values of r starting from 
the initial condition xq -0,Vo- 0. As t approaches t* , the 
probability density near x =; -0.4 increases monotonically. 
For T > T* the trend is reversed. Furthermore, it can be 
checked using Eq. (20) that for e = -0.4, I{x) attains its 
maximum around x -0.4. It is worth stressing again that 
even in this case, due to the extremely long time scales 
for small bias u, an observation up to would not show 
the reaching of a common steady state for different ini¬ 
tial conditions of the oscillator and may even lead to an 
apparent divergence of the {/(t)> traces. 
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Figure 10 Reduced probability density d^{x-,r) for u = 0.2, 
£ = -0.4 and different values of t, given the initial condition 
xo = 0,vo = 0- T = 150 (red); t = 1000 (green); r = 5.5- 10^ 
(blue); T = 10^ (purple); t = 4-10^ (cyan). Other parameters: 
0) - 10~^, 7= 0.08. 


3.4 Discussion 

3.4.1 Validity of the model and of the strong coupling 
regime 

The main assumption of this paper is the adiabatic con¬ 
dition 0 ) « 7 , which allows to neglect the dot transient 
dynamics in Eqns. (8,10,9) and ensuing Eqns. (12,14,13) 
which only retain the oscillator time dependence x{t). In¬ 
deed, this approximation remains valid as long as the adi¬ 
abatic condition is fulfilled and one explores time scales 
slower ( 1 )^^. Increasing oj ~ j breaks the validity of the 
above approximation [36]. In this case, it can be expected 
that qualitative modifications can occur with respect to 
the short-time behavior described in this paper. However, 
such modifications should not spoil the uniqueness of the 
steady state. This interesting topic is beyond the scope of 
the present paper and may be the subject of future inves¬ 
tigations. 

Additionally, in order to observe bistability and tristability, 
the strong coupling condition 7 < 1 , w « 1 must be ful¬ 
filled [32], as mentioned in Sec. 2. 

Quantum fluctuations beyond the damping and noise 
terms in the semiclassical Langevin equation may in 
principle be present [21,65]. Such fluctuations may as 


well, influcence the transient dynamics of the system 
affecting the hopping probability of the oscillator state. 
However, they have been shown to be relevant when 
u<(D [29,32,37]. Thus, all results presented here are well 
within the domain of the semiclassical approximation. 

Thermal flutuations induced by an extrinsic thermal 
bath as well as intrinsic mechanical damping mecha¬ 
nisms have been neglected in this paper. They can be 
modeled by an additional (constant) damping term in 
A{x] — A{x) + rj and diffusion term D{x] — Dix) + rjTri 
where = kg Tg/ZEp and Tg is the effective bath temper¬ 
ature [32]. Both the enhancement of damping and hopp- 
ping are expected to quantitatively affect the time scales 
discussed in the paper. In particular, thermal fluctuations 
are expected to lead, in the regime of strong coupling to 
this additional effective thermal bath, to an enhancement 
of the hopping rate and to a shortening of the dwell time 
Trf. Physically, we can expect that our results are not qual¬ 
itatively affected if t] < Amax and t? Tjj < Hmax where we 
can estimate Amax ~ ojlj^ and Hmaxaw/y For the param¬ 
eters employed in this paper, one gets Amax 10 ^^ and 
Fornax ~ 10 

Let us now relate these conditions to experimental pa¬ 
rameters. Supposing that for ri = aj/j^ = A^ax the damp¬ 
ing is essentially dominated by the extrinsic bath and 
still assuming a weak renormalization of the system fre¬ 
quency due to damping, one can identify the quality fa.c- 
tor Q =2 . Thus, the condition r] < Amax is equivalent to 

Q » Amax- For the case considered in this paper, Q > 10 is 
therefore required. As will be seen, systems such as sus¬ 
pended CNTs can easily exceed this limitation. 

Assuming the worst-case scenario rj = w/j^, the condi¬ 
tion T]Tp < Hmax can be cast into a bound for Tp which 
implies kg Th<To. Assuming a typical To =2 10 GHz one 
has Th<0.l K. 


3.4.2 Observability of the proposed effect 

The fact that the oscillator properties ({x(t)>) are directly 
related to those of the dot ({/(t))) make our results par¬ 
ticularly appealing, since the latter should in principle be 
detectable much more easily in a transport experiment. 
As already discussed above, in order to detect the dwell 
and relaxation time scales and Xr, one needs to con¬ 
sider a strong coupling between the electrons and the 
adiabatic oscillator. In addition to the requirements about 
the quality factor and extrinsic thermal bath tempera¬ 
tures, further restrictions apply on the bare system oscil¬ 
lator frequency, tunneling rate and the electron-vibration 
coupling force. In the following we will review three of 
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the most appealing candidates to observe the transient 
dynamics, namely nano-beams [28], suspended carbon 
nanotubes [5] and molecular systems [66,67] 
Nano-mechanical cantilevers, such as those that may be 
created in Si structures [28] are characterized by typical 
bare vibrational frequencies as low as 1-r 10 MHz, which 
would yield 1 -r 10 ms and Xr ~ 10 -r 100 ms at small 
bias. Similar performances may be obtained with metal 
nano-beams coated with a semiconducting piezoresistive 
layer [3]. With such low vibrational frequencies, satisfy¬ 
ing the adiabatic condition poses no particular challenge. 
The most relevant coupling mechanism between elec¬ 
trons and the oscillations is due to electrostatic gating. 
The capacitance of the beam with respect to an external 
gate fluctuates as the flexural mode is excited, leading to a 
bilinear coupling between the excess charge on the nano¬ 
beam and the amplitude of the flexural mode. Estimates 
show that such coupling is not very large [28], although it 
can be in principle tuned by a suitable a gate voltage. 
Another possible candidate is a suspended carbon nan¬ 
otube [5] considering the lowest flexural mode. Such a 
system can easily satisfy the adiabatic condition and ex¬ 
hibit strong electron-vibron coupling. In view of the great 
interest and relevance of this system, we provide here ac¬ 
tual estimates based on state-of-the-art experimental se¬ 
tups. We begin discussing the coupling between electrons 
and vibrations, which occurs in two different ways. The 
mtnnsic electron-phonon mechanism [68] would lead to 
a coupling term quadratic in the displacement operator 
X. This coupling, however, has been shown to be very 
weak [68] and unlikely to give rise to bistability. However, 
also in this case an external biased gate can be employed 
to couple the flexural mode and electrons. In this case, the 
coupling is linear in X and has the form of Eq. (3) with the 
electron-phonon coupling force given by [69] 


(Co -I- Cg)ln [d/r] d ’ 


(32) 


where 
^ _ ZnepL 
In (d/r) 


(33) 


is the CNT-gate capacitance, with cp the vacuum permit¬ 
tivity, L the CNT length, d the distance between the CNT 
and the gate and r the CNT radius. Finally, Co is the ca¬ 
pacitance of the CNT with respect to all other gates, and 
Vg the voltage present on the gate. 

We will consider a CNT with typical parameters L ^ 
1.8 pm, r =! 1.5 nm, d =: 400 nm and a typical gate volt¬ 
age Vg 1 V [70]. We also assume a typical value of 
Co 10^^^ F [72]. One obtains Ho 200 MHz =; 0.8 peV 
where a stiffening contribution due to the additional 


strain induced by the gate has been accounted for [71], 
Cg ~ 2 ■ 10“^^ F and thus Ep = A^/ZmD.^ ^ 15 peV, with 
m ^ 10“^*’ Kg. This implies w =; 0.02 « 1. Setting y = 
lOta = 0.2 < 1 implies a typical tunneling rate To ~ 6 GHz 
corresponding to an average current Ip~l nA, well within 
the range of experiments [70]. With the above estimates 
we get Trf 0.1 ms and ~ 1ms. Extremely high quality 
factors Q > 10^ have been reported for the flexural mode 
of suspended CNTs [73], constituting a favorable condi¬ 
tion to observe the predicted effects. Other vibrational 
modes or tensile-stressed CNTs [10] have frequencies at 
least two orders of magnitude larger [6], making the de¬ 
tection of the above times more problematic. 

Finally, also molecular systems exhibit a strong coupling 
between electrons and vibrations [32]. However, the cou¬ 
pling to the leads is generally very weak, which in addi¬ 
tion to the usually large oscillation frequencies makes it 
extremely difficult to achieve the adiabatic regime [32,33] 
considered in this paper. 


4 Conclusions 


In this paper we have investigated the time evolution to¬ 
wards the steady state of a NEMS in the adiabatic regime. 
The Eangevin equation for the oscillator and the corre¬ 
sponding Fokker-Planck equation for the probability den¬ 
sity in phase space are numerically solved for arbitrary 
time and initial conditions in several parameter settings 
of the system, corresponding to bi-stability (low bias), 
tri-stability (intermediate bias and resonance) or mono- 
stability (large bias). It is shown that in all situations a 
unique steady state is reached. The approach to the steady 
state is however different according to the considered pa¬ 
rameter setting. For low bias, the dynamics of the system 
is found to be essentially frozen in a sort of blockaded 
regime up to rather large time scales, before relaxing to¬ 
wards the steady state. Two time scales emerge, corre¬ 
sponding to the dwell time into quasi-stable minima of 
the effective potential well of the oscillator and to the 
relaxation time towards the steady situation. For inter¬ 
mediate bias, the dynamics is described by either one 
or two time scales, depending on the initial condition of 
the system. For large bias, the system relaxes towards the 
steady state with a single time scale. Dwell and relaxation 
time scales have been identified numerically by analyzing 
the statistical properties of the solutions of the Eangevin 
equation and the spectral properties of the Fokker-Planck 
operator. Our findings suggest that no hysteretic behav¬ 
ior is shown by this system at steady state, but that an 
observation of the system dynamics stopping at short 
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time scales might fail to exhibit the tendency towards a 
unique steady state, especially at small bias and/or off- 
resonance [55]. The predicted time scales are within the 
range of current experiments. 
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